5 - Factor analysis

Author

CDN team

Published

March 1, 2024

Introduction

In this notebooks we are going to carry out Factor Analysis analysis using RcppML. You can see the GitHub repository here. This implementation of NMF is extremely fast and enables us to use large dataset since it works well with sparse matrices. The author is currently implementing it within the Singlet package to work nicely with single cell and soon it will be available!

here and the differential communication analysis vignette here. The goal of the notebook is to identify which communication pathways are altered between flu and covid infection!

Some associated literature which is a must read are:

Key Takeaways

  • NMF identifies sets of genes “metagenes” representing the main characteristics of the data.

  • Choosing the rank (K) of the matrix is an important step since it will determine how many “metagenes” are present in our dataset. We can use cross validation to find the best K.

  • Whe can visualize how important each genes is for each factor (aka - metagene) and how important is each factor for each cell. This way we can determine which metagenes are explaining each cell’s transcriptome.

  • When looking at metagenes learned across multiple patients and conditions it is imperative to check that the signal is not being provided just by one replicate of our experiment. Not checking this could lead to incorrect interpretation of the data!

Library

options(future.globals.maxSize= 891289600)
### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
    install.packages("Seurat")

if (!requireNamespace("dplyr", quietly = TRUE))
    install.packages("dplyr")

if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!requireNamespace("RcppML", quietly = TRUE)) {
    # BiocManager::install("fgsea", update = FALSE)
    # BiocManager::install("limma", update = FALSE)
    # devtools::install_github("zdebruine/singlet", upgrade = FALSE)
    devtools::install_github("zdebruine/RcppML")
    devtools::install_github("zdebruine/RcppSparse")
}

if (!requireNamespace("sparseMatrixStats", quietly = TRUE)) {
    install.packages("sparseMatrixStats")
}

if (!requireNamespace("inflection", quietly = TRUE)) {
    install.packages("inflection")
}

if (!requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
    BiocManager::install("org.Hs.eg.db")
}

if (!requireNamespace("clusterProfiler", quietly = TRUE)) {
    BiocManager::install("clusterProfiler")

}


### Load all the necessary libraries
library(Seurat)
library(dplyr)
library(RcppML)
library(RcppSparse)
library(ggplot2)
library(sparseMatrixStats)
library(inflection)
library(glue)
library(tidyr)
library(org.Hs.eg.db)
library(clusterProfiler)

set.seed(687)

Load data

We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.

# Download the data in data/ directory
# download.file(
#     url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
#     destfile = "../data/workshop-data.rds",
#     method = "wget",
#     extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds

se <- readRDS("../data/workshop-data.rds")

Analysis

Convert ENSEMBL IDs to Gene Symbols

Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:

head(rownames(se))
[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
[5] "ENSG00000000460" "ENSG00000000938"

Convert to gene symbols

gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")

symbol_id <- data.frame(index = rownames(se)) %>%
    left_join(gene_df) %>%
    pull(feature_name)

# re-create seurat object
mtx <- se@assays$RNA@data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
# symbol_id <- mapIds(
#     org.Hs.eg.db,
#     keys = rownames(se), 
#     column = "SYMBOL",
#     keytype = "ENSEMBL",
#     multiVals = "first")
# 
# # df <- data.frame(symbol = symbol_id, ensembl = names(symbol_id))
# all(rownames(se) == names(symbol_id))
# 
# # re-create seurat object
# mtx <- se@assays$RNA@data
# rownames(mtx) <- symbol_id
# 
# # Remove NAs
# sum(is.na(rownames(mtx)))
# dim(mtx)
# mtx <- mtx[!is.na(rownames(mtx)), ]
# dim(mtx)
# 
# rownames(mtx) <- make.unique(rownames(mtx))
# se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
# 
# sum(is.na(rownames(se)))

Factor Analysis

For the purpose of this vignette we’re going to focus on CD8 and NK cells this object to pull out interferon signalling program:

se <- se[, se$Celltype %in% c("CD8, non-EM-like", "CD8, EM-like", "NK cell")]

First we need to normalize the data

se <- NormalizeData(se, verbose = FALSE) %>%
    FindVariableFeatures(nfeatures = 3000, verbose = FALSE) %>%
    ScaleData(verbose = FALSE) %>%
    RunPCA(verbose = FALSE)

VlnPlot(
    se,
    features = c("CD3D", "CD3E", "CD8B", "NKG7", "FCGR3A"),
    group.by = "Celltype",
    pt.size = 0,
    split.by = "disease") + theme(legend.position = "bottom")

Extract the normalized expression matrix and remove genes that are all 0s

library(Matrix)
mtx <- se@assays$RNA@layers$data
colnames(mtx) <- colnames(se)
rownames(mtx) <- rownames(se)

# Remove genes that are all 0s
mtx <- mtx[sparseMatrixStats::rowSums2(mtx) > 0, ]

Cross-validation for rank determination

The first step is to determine the rank of our matrix - by this we mean, which is the optimal number of factors we need to decompose our matrix. Here we run cross-validation 3 times across a range of different ranks (k).

start <- Sys.time()
cv_data <- crossValidate(mtx, k = seq(1, 36, 5), reps = 2, verbose = TRUE, tol = 1e-04)
Sys.time() - start
Time difference of 4.203059 mins

Now let’s visualize the crossvalidation results

ggplot(cv_data, aes(x = k, y = value, color = rep, group = rep)) +
    geom_point() +
    geom_line() +
    theme_minimal()

We can see how it starts to plateau at ~k=20 so we will use that in this script.

Run NMF

  • k: Number of factors we want to identify

  • tol: Stands for tolerance and is how small we want the error to be, the smaller the number the better the decomposition but the longer its gonna take to converge. Values in the 10^-5 and 10^-6 return very good results.

  • L1: Introduces sparsity into the factors and loadings with the aim of removing noisy genes from factors and noisy cells contributing to factors. L1 normalization uses Lasso regularization to increase sparsity and remove the unimportant features.

set.seed(7)
nmf_ls <- RcppML::nmf(
    data = mtx,
    k = 20,
    tol = 1e-06,
    L1 = c(0.05, 0.1), # l1 for c(w, h)
    verbose = TRUE)

Explore the NMF results

W - contains which genes are relevant for each factor. H - contains how important each factor for each cell.

w <- nmf_ls@w
w[1:5, 1:5]
         nmf1         nmf2         nmf3         nmf4         nmf5
TSPAN6      0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
DPM1        0 9.446519e-05 4.755987e-05 1.160628e-04 7.280245e-05
SCYL3       0 0.000000e+00 0.000000e+00 5.235962e-07 0.000000e+00
C1orf112    0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
FGR         0 5.544984e-05 4.062468e-04 5.537400e-05 2.077896e-04
h <- nmf_ls@h
h[1:5, 1:5]
     AAACCCAAGAATGTTG-12 AAACCCAAGCTACGTT-6 AAACCCAAGGCCTGCT-12
nmf1        5.462850e-05       5.187584e-05        3.464711e-05
nmf2        7.132754e-05       1.925005e-04        0.000000e+00
nmf3        1.210247e-04       2.170797e-05        0.000000e+00
nmf4        1.459895e-04       0.000000e+00        0.000000e+00
nmf5        5.263035e-05       0.000000e+00        1.276532e-04
     AAACCCAAGTGTAGAT-5 AAACCCACACAATGCT-5
nmf1       6.166142e-05       5.885350e-05
nmf2       9.268590e-05       1.126141e-04
nmf3       0.000000e+00       0.000000e+00
nmf4       4.154653e-06       2.351866e-05
nmf5       7.246126e-05       1.722137e-05

Add NMF to Seurat object

se[["FA"]] <- Seurat::CreateDimReducObject(
  embeddings = t(nmf_ls$h),
  loadings = nmf_ls$w,
  assay = "RNA",
  key = "Factor_")

se[["FA"]]
A dimensional reduction object with key Factor_ 
 Number of dimensions: 20 
 Number of cells: 19262 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

Examine and visualize NMF loadings a few different ways

print(se[["FA"]], dims = 1:20, nfeatures = 5)
Factor_ 1 
Positive:  MALAT1, MT-CO3, MT-CO2, RPL10, EEF1A1 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 2 
Positive:  GNLY, CD52, GZMH, FGFBP2, NKG7 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 3 
Positive:  GNLY, GZMB, NKG7, GZMH, IL32 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 4 
Positive:  GNLY, FGFBP2, NKG7, TYROBP, CCL5 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 5 
Positive:  NFKBIA, GNLY, CD69, CCL4, IER2 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 6 
Positive:  JUN, JUNB, CXCR4, DUSP1, TSC22D3 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 7 
Positive:  S100A9, S100A8, PRF1, NKG7, GZMA 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 8 
Positive:  CMC1, CD3D, NKG7, HLA-DRB1, CD74 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 9 
Positive:  ACTB, LGALS1, CALR, HSP90AB1, PFN1 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 10 
Positive:  MALAT1, MT-CO2, MT-CO3, MT-CO1, MT-ATP6 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 11 
Positive:  CREM, METRNL, CXCR4, DUSP2, PIK3R1 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 12 
Positive:  FCER1G, SPON2, MYOM2, FGFBP2, IGFBP7 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 13 
Positive:  CREM, ZNF331, AREG, CEMIP2, DUSP2 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 14 
Positive:  CD74, ACTG1, HLA-DRA, COTL1, LIMD2 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 15 
Positive:  GNLY, CTSW, SELL, GZMK, CD7 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 16 
Positive:  ISG15, MTRNR2L12, XAF1, JUN, IFI6 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 17 
Positive:  LTB, IL7R, CD52, VIM, IL32 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 18 
Positive:  RPS26, METRNL, HLA-DQB1, HLA-DRB1, ID2 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 19 
Positive:  GNLY, GZMB, PRF1, GZMH, FLNA 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 
Factor_ 20 
Positive:  KLRD1, TXNIP, ZFP36L2, KLF2, RPS4Y1 
Negative:  TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7 

Violin plots for the loadings

VlnPlot(
  se,
  features = glue::glue("Factor_{1:ncol(nmf_ls$w)}"),
  ncol = 4,
  group.by = "Celltype")

We can also visualize this by cell type, condition and patient!

disease_pal <- c("#41AE76", "#225EA8", "#E31A1C")
names(disease_pal) <- c("normal", "influenza", "COVID-19")
# flu <- RColorBrewer::brewer.pal(12, name = "YlGnBu")
# normal <- RColorBrewer::brewer.pal(12, name = "BuGn")
# covid <- RColorBrewer::brewer.pal(12, name = "YlOrRd")
donor_pal <- c(
    "#66C2A4", "#41AE76", "#238B45", "#006D2C",
    "#41B6C4", "#1D91C0", "#225EA8", "#253494", "#081D58",
    "#FFFFCC", "#FFEDA0", "#FED976", "#FEB24C", "#FD8D3C",
    "#FC4E2A", "#E31A1C", "#BD0026", "#800026")

names(donor_pal) <- c(
    "Normal 1", "Normal 2", "Normal 3", "Normal 4",
    "Flu 1", "Flu 2", "Flu 3", "Flu 4", "Flu 5",
    "nCoV 1", "nCoV 2", "nCoV 3_4", "nCoV 5",
    "nCoV 6", "nCoV 7_8", "nCoV 9_10", "nCoV 11"  
)

# Preprocess dataset
dd <- bind_cols(se@meta.data, se@reductions$FA@cell.embeddings) %>%
        tidyr::pivot_longer(
            cols = glue::glue("Factor_{1:ncol(nmf_ls$w)}"),
            names_to = "k",
            values_to = "loading") %>%
        group_by(Celltype, disease, donor_id, k) %>%
        summarise(median_loading = median(loading))

lapply(glue::glue("Factor_{1:ncol(nmf_ls$w)}"), function(i) {
    dd %>%
        filter(k == i) %>%
        ggplot(aes(x = disease, y = median_loading, fill = disease)) +
        geom_boxplot(outlier.shape = NA, alpha = 0.7) +
        geom_jitter(aes(color = donor_id), height = 0) +
        facet_wrap(~ Celltype) +
        scale_fill_manual(values = disease_pal) +
        scale_color_manual(values = donor_pal) +
        scale_x_discrete(limits = names(disease_pal)) +
        theme_classic() +
        labs(title = glue::glue("Factors by condition in {i}"))
})
[[1]]


[[2]]


[[3]]


[[4]]


[[5]]


[[6]]


[[7]]


[[8]]


[[9]]


[[10]]


[[11]]


[[12]]


[[13]]


[[14]]


[[15]]


[[16]]


[[17]]


[[18]]


[[19]]


[[20]]

Factors 4 & 6 seem interesting, lets explore them a bit more in depth!

w[1:5, c("nmf5", "nmf7")]
                 nmf5         nmf7
TSPAN6   0.000000e+00 0.0000000000
DPM1     7.280245e-05 0.0000000000
SCYL3    0.000000e+00 0.0000000000
C1orf112 0.000000e+00 0.0000000000
FGR      2.077896e-04 0.0005231307
# Extract top genes from factors 4 & 6 using Unit Invariant Knee
FA_genes <- lapply(c("nmf5", "nmf7"), function(i) {
    y_vec <- sort(w[, i], decreasing = TRUE)
    # Define inflecion point
    n <- uik(x = seq_len(length(y_vec)), y = y_vec)
    # Define top 10 genes
    df <- data.frame(gene = names(y_vec), value = y_vec, rank = 1:length(y_vec)) %>%
        mutate(lab = if_else(rank < 15, gene, NA_character_))
    print(ggplot(df, aes(x = rank, y = y_vec)) +
        geom_point() +
        geom_vline(xintercept = n) +
        ggrepel::geom_text_repel(aes(label = lab)) +
        labs(title = glue("Factor-{i}")) +
        theme_classic())
    y_vec[1:n]
})

names(FA_genes) <- c("nmf5", "nmf7")

dd <- lapply(names(FA_genes), function(i){
    data.frame(
        value = FA_genes[[i]],
        gene = names(FA_genes[[i]]),
        factor = i)
}) %>%
    bind_rows() %>%
    pivot_wider(names_from = factor, values_from = value, values_fill = 0)

DT::datatable(dd)

By looking at these genes we can get a sense of which are the major processes captures by each factor.

  • NMF-5 seems to be capturing cytotoxic and activation signals since it contains genes encoding for chemokines (CCL4, CCL3), granzymes (GZMB, GNLY, NKG7) and other immune related processes (NFKBIA, CD69…).

  • NMF-7 contains a wide array of myeloid cell markers, such as S100A8, S100A9, LYZ and genes encoding for MHC-II (CD74, HLA-DRB1, HLA-DRA) as well as cytotoxic genes (GZMA, GZMH, NKG7…). The most straight forward explanation is that these could be doublets or have a high amount of ambient RNA in the soup.

We can take a look at how these genes are expressed across the patients since like Flu-1 patient has an etremely high score and we want to make sure it is not driving this by itself.

# We'll take a look at the top 25 genes
names(head(FA_genes[["nmf7"]], 25))
 [1] "S100A9"  "S100A8"  "PRF1"    "NKG7"    "GZMA"    "IL32"    "CD74"   
 [8] "GZMB"    "LYZ"     "GZMH"    "HLA-B"   "B2M"     "HBB"     "TXNIP"  
[15] "S100A11" "MT-CO1"  "PSME2"   "PSMB9"   "PLAAT4"  "RPS4Y1"  "CST7"   
[22] "TRAC"    "FGFBP2"  "ACTB"    "HLA-A"  
# Scale their expression
se <- ScaleData(se, features = names(head(FA_genes[["nmf7"]], 25)))

# Check scale.data
mtx <- se@assays$RNA@layers$scale.data
rownames(mtx) <- names(head(FA_genes[["nmf7"]], 25))
colnames(mtx) <- colnames(se)

# Choose color palette
palette_length <- 100
my_color <- colorRampPalette(c("Darkblue", "white","red"))(palette_length)

# set color breaks
my_breaks <- c(seq(quantile(mtx, .01), 0, length.out=ceiling(palette_length/2) + 1),
               seq(0.05, quantile(mtx, .99), length.out=floor(palette_length/2)))
# Note that we are setting the min and max of the scale to .01 and .99 so we exclude extreme values from dampening our signal
# We can see how the max value of the matrix is +10 but q99 is 2
Hmisc::describe(as.vector(mtx))
as.vector(mtx) 
         n    missing   distinct       Info       Mean        Gmd        .05 
    481550          0     287115          1 -2.226e-06        1.1   -1.59884 
       .10        .25        .50        .75        .90        .95 
  -1.31044   -0.59795    0.02677    0.69204    1.16664    1.48214 

lowest : -12.5361 -8.98419 -8.74596 -8.11004 -7.98887
highest: 9.87275  9.88233  9.96698  9.97462  10      
pheatmap::pheatmap(
    mtx,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = FALSE,
    treeheight_col = 0,
    annotation_col = se@meta.data[, c("disease", "donor_id", "Celltype")], 
    annotation_colors = list("disease" = disease_pal, "donor_id" = donor_pal),
    color = my_color,
    breaks = my_breaks)

Carry out GSEA

# read GSEA markers
gsea_ls <- lapply(FA_genes, function(i) {
    
    # http://yulab-smu.top/clusterProfiler-book/chapter5.html#go-gene-set-enrichment-analysis
    gsea_results <- clusterProfiler::gseGO(
        geneList = i,
        ont = "BP",
        OrgDb = org.Hs.eg.db,
        keyType = "SYMBOL",
        minGSSize = 10,
        maxGSSize = 300,
        pvalueCutoff = 0.1,
        pAdjustMethod = "BH",
        seed = TRUE)
    
    gsea_results
})

names(gsea_ls) <- names(FA_genes)

Visualize top enriched gene sets per cluster

lapply(names(gsea_ls), function(i) {
    # Extract gsea
    gsea <- gsea_ls[[i]]
    gsea <- clusterProfiler::simplify(gsea, cutoff = 0.7)
    gsea@result <- gsea@result %>%
        dplyr::arrange(dplyr::desc(NES)) %>%
        dplyr::filter(p.adjust < 0.1)
    
    tmp_plt <- gsea@result %>%
        dplyr::top_n(n = 20, wt = enrichmentScore) %>%
        ggplot2::ggplot(.,
                        ggplot2::aes(
                            x = NES,
                            y = forcats::fct_reorder(Description, NES),
                            size = setSize,
                            color = p.adjust)) +
        ggplot2::geom_point() +
        # ggplot2::geom_vline(
        #     xintercept = 0,
        #     linetype = "dashed",
        #     color = "red") +
        scale_color_viridis_c(option = "plasma") +
        ggplot2::theme_minimal() +
        ggplot2::labs(title = i)
    

  }) %>% patchwork::wrap_plots(ncol = 2)

And lastly we will save GSEA to a spreadsheet so we can take a look:

gsea_xlsx <- lapply(names(gsea_ls), function(i) {
  print(i)
  # Extract gsea
  gsea <- gsea_ls[[i]]
  gsea <- clusterProfiler::simplify(gsea, cutoff = 0.7)
  gsea@result <- gsea@result %>%
    dplyr::arrange(dplyr::desc(NES)) %>%
    dplyr::filter(p.adjust < 0.1)
  })
[1] "nmf5"
[1] "nmf7"
names(gsea_xlsx) <- names(gsea_ls)
openxlsx::write.xlsx(gsea_xlsx, file = "../data/GSEA.xlsx")

Session Info

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] Matrix_1.6-5             clusterProfiler_4.8.2    org.Hs.eg.db_3.17.0     
 [4] AnnotationDbi_1.62.2     IRanges_2.34.1           S4Vectors_0.38.1        
 [7] Biobase_2.60.0           BiocGenerics_0.46.0      tidyr_1.3.1             
[10] glue_1.7.0               inflection_1.3.6         sparseMatrixStats_1.12.2
[13] MatrixGenerics_1.12.3    matrixStats_1.2.0        ggplot2_3.5.0           
[16] RcppSparse_1.0           RcppML_0.5.6             dplyr_1.1.4             
[19] Seurat_5.0.1             SeuratObject_5.0.1       sp_2.1-3                

loaded via a namespace (and not attached):
  [1] fs_1.6.3                spatstat.sparse_3.0-3   bitops_1.0-7           
  [4] enrichplot_1.20.0       devtools_2.4.5          HDO.db_0.99.1          
  [7] httr_1.4.7              RColorBrewer_1.1-3      backports_1.4.1        
 [10] profvis_0.3.8           tools_4.3.1             sctransform_0.4.1      
 [13] DT_0.29                 utf8_1.2.4              R6_2.5.1               
 [16] lazyeval_0.2.2          uwot_0.1.16             urlchecker_1.0.1       
 [19] withr_3.0.0             prettyunits_1.1.1       gridExtra_2.3          
 [22] progressr_0.14.0        cli_3.6.2               spatstat.explore_3.2-6 
 [25] fastDummies_1.7.3       scatterpie_0.2.1        sass_0.4.8             
 [28] labeling_0.4.3          spatstat.data_3.0-4     readr_2.1.4            
 [31] ggridges_0.5.6          pbapply_1.7-2           yulab.utils_0.1.4      
 [34] foreign_0.8-84          gson_0.1.0              DOSE_3.26.2            
 [37] parallelly_1.37.0       sessioninfo_1.2.2       rstudioapi_0.15.0      
 [40] RSQLite_2.3.1           generics_0.1.3          gridGraphics_0.5-1     
 [43] crosstalk_1.2.1         ica_1.0-3               spatstat.random_3.2-2  
 [46] vroom_1.6.3             zip_2.3.0               GO.db_3.17.0           
 [49] ggbeeswarm_0.7.2        fansi_1.0.6             abind_1.4-5            
 [52] lifecycle_1.0.4         yaml_2.3.8              qvalue_2.32.0          
 [55] Rtsne_0.17              grid_4.3.1              blob_1.2.4             
 [58] promises_1.2.1          crayon_1.5.2            miniUI_0.1.1.1         
 [61] lattice_0.21-8          cowplot_1.1.3           KEGGREST_1.40.0        
 [64] pillar_1.9.0            knitr_1.45              fgsea_1.26.0           
 [67] future.apply_1.11.1     codetools_0.2-19        fastmatch_1.1-4        
 [70] leiden_0.4.3.1          downloader_0.4          ggfun_0.1.4            
 [73] data.table_1.15.0       remotes_2.4.2.1         vctrs_0.6.5            
 [76] png_0.1-8               treeio_1.24.3           spam_2.10-0            
 [79] gtable_0.3.4            cachem_1.0.8            openxlsx_4.2.5.2       
 [82] xfun_0.42               mime_0.12               tidygraph_1.2.3        
 [85] survival_3.5-7          pheatmap_1.0.12         ellipsis_0.3.2         
 [88] fitdistrplus_1.1-11     ROCR_1.0-11             nlme_3.1-163           
 [91] ggtree_3.8.2            usethis_2.2.2           bit64_4.0.5            
 [94] RcppAnnoy_0.0.22        GenomeInfoDb_1.36.3     bslib_0.6.1            
 [97] irlba_2.3.5.1           rpart_4.1.19            vipor_0.4.5            
[100] KernSmooth_2.23-22      Hmisc_5.1-0             colorspace_2.1-0       
[103] DBI_1.1.3               nnet_7.3-19             ggrastr_1.0.2          
[106] tidyselect_1.2.0        processx_3.8.2          bit_4.0.5              
[109] compiler_4.3.1          htmlTable_2.4.1         plotly_4.10.4          
[112] shadowtext_0.1.3        checkmate_2.2.0         scales_1.3.0           
[115] lmtest_0.9-40           callr_3.7.3             stringr_1.5.1          
[118] digest_0.6.34           goftest_1.2-3           spatstat.utils_3.0-4   
[121] rmarkdown_2.25          XVector_0.40.0          base64enc_0.1-3        
[124] htmltools_0.5.7         pkgconfig_2.0.3         fastmap_1.1.1          
[127] rlang_1.1.3             htmlwidgets_1.6.4       shiny_1.8.0            
[130] jquerylib_0.1.4         farver_2.1.1            zoo_1.8-12             
[133] jsonlite_1.8.8          BiocParallel_1.34.2     GOSemSim_2.26.1        
[136] RCurl_1.98-1.12         magrittr_2.0.3          Formula_1.2-5          
[139] GenomeInfoDbData_1.2.10 ggplotify_0.1.2         dotCall64_1.1-1        
[142] patchwork_1.2.0         munsell_0.5.0           Rcpp_1.0.12            
[145] ape_5.7-1               viridis_0.6.4           reticulate_1.35.0.9000 
[148] stringi_1.8.3           ggraph_2.1.0            zlibbioc_1.46.0        
[151] MASS_7.3-60             plyr_1.8.9              pkgbuild_1.4.2         
[154] parallel_4.3.1          listenv_0.9.1           ggrepel_0.9.5          
[157] forcats_1.0.0           deldir_2.0-2            Biostrings_2.68.1      
[160] graphlayouts_1.0.0      splines_4.3.1           tensor_1.5             
[163] hms_1.1.3               ps_1.7.5                igraph_2.0.2           
[166] spatstat.geom_3.2-8     RcppHNSW_0.6.0          reshape2_1.4.4         
[169] pkgload_1.3.2.1         evaluate_0.23           BiocManager_1.30.22    
[172] tzdb_0.4.0              tweenr_2.0.2            httpuv_1.6.14          
[175] RANN_2.6.1              purrr_1.0.2             polyclip_1.10-6        
[178] future_1.33.1           scattermore_1.2         ggforce_0.4.1          
[181] xtable_1.8-4            RSpectra_0.16-1         tidytree_0.4.6         
[184] later_1.3.2             viridisLite_0.4.2       tibble_3.2.1           
[187] aplot_0.2.2             memoise_2.0.1           beeswarm_0.4.0         
[190] cluster_2.1.4           globals_0.16.2